ELECTROSTATIC FORCES ON CHARGED SURFACES OF 
BILAYER LIPID MEMBRANES * 



MICHAEL MIKUCKF AND Y. C. ZHOU* 

Abstract. Simulating protein-membrane interactions is an important and dynamic area of 
research. A proper definition of electrostatic forces on membrane surfaces is necessary for developing 
electromechanical models of protein-membrane interactions. Here we modeled the bilayer membrane 
as a continuum with general continuous distributions of lipids charges on membrane surfaces. A 
new electrostatic potential energy functional was then defined for this solvatcd protein-membrane 
system. We investigated the geometrical transformation properties of the membrane surfaces under a 
smooth velocity field. These properties allows us to apply the Hadamard-Zolcsio structure theorem, 
and the electrostatic forces on membrane surfaces can be computed as the shape derivative of the 
' electrostatic energy functional. 

a ; 

Key words, lipid bilayer membrane; electromechanics; dielectric interface; surface charge; 
■ dielectric boundary force; implicit-solvent; Poisson-Boltzmann; shape derivative 



o 



PL, 



43 



o 



o 



- 1—1 

X 



AMS subject classifications. 35J, 35Q, 49S, 82B, 82D, 92C 



1. Introduction. This paper concerns the mathematically rigorous and physi- 
cally justifiable definition of electrostatic forces on the surfaces of lipid bilayer mem- 
branes within the framework of implicit solvent and continuum models of charged 
lipids. Lipid bilayer membranes, as the boundaries of cells and many cell organelles, 
control the exchange of ions, nutrient particles and metabolic products between the 
enclosed structures and the surrounding aqueous environment. Bilayer membranes 
function by stretching, bending, merging and separating to control gate specific chan- 
nels or to wrap/unwrap the particles. These deformations are precisely regulated 
by various proteins and other macromolecules, each with its own specific functions. 
Among all types of intermolecular interactions that drive the membrane deformation, 
the electrostatic interaction is ubiquitous because proteins and lipid bilayer mem- 
branes are always charged under physiological conditions. It is indeed one of the 
most important interactions if the membrane dynamics involve the lateral diffusion 
(•r-j ■ of charged lipids and electrostatic association of proteins. Computation of the elec- 

trostatic forces on the membrane is therefore of critical importance for quantitative 
study of membrane dynamics and related cellular activities. Interested readers are 
referred to [13, 21, 18, 3, 17] for more thorough discussions of the membrane dynamics 
and protein-membrane electrostatic interactions. 

Computation of electrostatic forces on bilayer membrane by summing up pair- 
wise Coulombic interactions has to consider full molecular details of the system and 
involves all particles in the computation. Thus, it is challenging to use this method 
to study membrane dynamics at biologically relevant spatial and time scales. Implicit 
solvent models have been introduced so that the averaged behavior of highly dynamic 
solvent molecules can be described as a structure- less continuum [9] . This simplifica- 
tion greatly reduces the degrees of freedom in simulations. A rigorous definition of the 
electrostatic force on biomolecules within the framework of implicit solvent models 
has been the subject of study for decades by computational biologists and applied 
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mathematicians [19, 11, 23, 16]. Recently, Li et al. obtained an elegant derivation 
of electrostatic forces on molecular surfaces by computing the shape derivative of 
the electrostatic potential energy [16]. That work assumes that the dielectric inter- 
face is uncharged. This assumption does not hold for bilayer membranes if they are 
modeled with continuous distribution of lipids on surfaces. Such a continuum mem- 
brane model is necessary for simulating membrane dynamics over a region larger than 
100A in space and longer than a microsecond in time [20, 4, 14, 24]. The variation 
of charge densities on dielectric interfaces has to be considered in deriving the elec- 
trostatic forces, and the model in [16] has to be modified before it can be applied 
to protein-membrane interactions. In this work, we propose an electrostatic energy 
functional for the Poisson-Boltzmann equation with a general variable charge density 
on bilayer membrane surfaces. Under a special condition this surface charge density 
follows a constraint Boltzmann distribution. We find that the time derivative of the 
Jacobian for surface transformation equals the surface divergence of the velocity field 
at membrane surfaces. Thanks to this essential geometric property, the Hadamard- 
Zolesio structure theorem holds true for our model and the shape derivative approach 
can be applied consequently to compute the electrostatic forces on bilayer membrane 
surfaces. 

The rest of the article is organized as follows. In Section 2 we introduce the 
governing equation and free-energy functional for the protein-membrane electrostatic 
interactions. Section 3 is devoted to the computation of the shape derivative of 
the electrostatic energy functional, with an emphasis on the treatment of membrane 
surface charge distribution and surface transformation. Finally, we shall show in 
Section 4 that the electrostatic force obtained through shape derivative approach 
matches the divergence of the jump of Maxwell stress tensor (MST). 

2. Motivating problem and mathematical model. The membrane defor- 
mation will be modeled under the following conditions. Define Vt C M 3 as the region 
containing the entire membrane-protein system. Let Q p denote the volume of the 
protein and fl m be the volume of the membrane. The solvent will occupy fi s which 
includes both the surrounding environment and the volume enclosed by the mem- 
brane. The boundary separating the protein Q p and the solvent fl s will be given by 
the manifold T p C R 2 . The membrane f2 m has two boundaries to the solvent. The 
interior boundary is called the cytosolic face and is denoted F c C M. 2 . The exte- 
rior boundary is called the endoplasmic face and is denoted F e C M 2 . The exterior 
boundary of the containment domain is called <9f2. The atomic detail of the lipids 
that compose the membrane will be neglected for the continuum model. The unit 
outward normal to any boundary F will be denoted as n. A cross section of the 
protein- membrane model is illustrated in Figure 2.1. 

The total potential energy of bilayer membrane is expressed in two components. 
First, the classical mechanical bending energy depends only on the curvature of the 
bilayer membrane, attributed to the work of Canham [2], Helfrich [12] and Evans [8]. 
Their calculation for the bending energy is given by 



where J£p and Jtc are the bending modulus and Gaussian saddle-splay modulus, 
respectively, H is the mean curvature, Co is the spontaneous curvature, K is the 
Gaussian curvature, a is the determinant of the covariant metric tensor, and s± and 
S2 are the intrinsic curvilinear coordinates of the membrane [10]. In addition to the 
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Fig. 2.1. Mathematical description of protein-membrane system. The containment domain is 
denoted by f!. The volume of the lipid bilayer is f! m , and the exdoplasmic (exterior) and cytosolic 
(interior) faces of the membrane are F e and T c , respectively. The volume enclosed by the membrane 
and the aqueous surrounding environment are both denoted f! s . The protein is Q p with a surface 
F p . The unit outward normal to any surface T is n. Lipids are drawn in the bottom left corner 
for illustration but their atomic details will be neglected in the model. Note that the protein may be 
located in the solvent region inside T c in some cases. 



bending energy, the membrane is under an external potential force induced by the 
protein. The potential energy from this interaction is added to the total potential 
energy of the system: 

n[r] = E[T)+G[T], (2.2) 

where E[T] is the bending energy (2.1) and G[T] is the electrostatic potential energy 
from the protein-membrane interaction. The equilibrium position of the bilayer mem- 
brane is determined by the surfaces which minimize both the bending energy and the 
electrostatic potential energy. The minimization of (2.2) gives rise to the bending 
curvature equation of T, 

S r Ii[T] = 5 r E[T] + 8 r G[T]. (2.3) 
The variation of the classical bending energy has been calculated in [10] to be 



S r E[T] = J 



JT c 2(2iJ - C Q ) 2 5Hy^ + J€ C -{2H - C ) 2 6y/a 



{1/VE) dS. (2.4) 



The main result of this article is the calculation of the external force induced by the 
protein, given by the variation of the electrostatic potential energy with respectto the 
boundary, i5rG[r]. 

2.1. Poisson-Boltzmann Equation with surface charge distributions. 

The electrostatic potential energy G[T] for the lipid membrane boundaries T e and 
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r c is related to the weak form of the nonlinear Poisson-Boltzmann equation, 



>]=0 
d4> s d(j/ 1 



dn 



dn 



£s dn £p dn 



in f2, 

on r c ,r e ,r p , 
on r c ,r e , 

on r„, 



(2.5) 



where e is the dielectric coefficient. Distinct dielectric permittivities are defined in 
Cl m , Cl s , and il p . The electrostatic potential is given by <fi. The subscripts of e and <\> 
denote the domain they are defined. The fixed charge density inside proteins is given 
by /, the charge of an individual lipid is given by qi, and the concentration of lipids 
on membrane surfaces is p[T]. Here, we assume there is only one species of diffusive 
charged lipids in the membrane, but the model can easily be generalized to multiple 
lipid species. The distribution of charged lipids on the membrane surfaces p\T] is 
governed by the surface electrodiffusion equation 



dt s 



{DV sP + DqipV s 



(2.6) 



where t is time, D is the diffusion coefficient, and V S ,V S - are the surface gradient 
and surface divergence operators, respectively. In this model, two distributions for the 
lipids p e = p[T e ] and p c = p[T c ] will be considered, ' corresponding to the solutions of 
(2.6) on the two membrane surfaces T e and T c , respectively. The function B describes 
the electrostatic energy due to the mobile ions and is defined by 



M 



-Pqj<P 



1), 



(2.7) 



where j3 = l/(kgT) is the inverse thermal energy, M is the number of ionic species 
in the solvent, and qj and c°° are the charge and bulk concentration of the j th ionic 
species, respectively [16]. 

One method of solving (2.6) utilizes a representation of the ion concentration p 
by a 'Slotboom variable' 



u = pe* 



(2.8) 



This change of variable transforms (2.6) to an equivalent, symmetric form. At the 
steady state, (2.6) becomes 



It is apparent that 



= V s • (L»e" w V s M). 



u = pe^ = c 



(2.9) 



(2.10) 



is a solution to (2.9) for some constant c. Then, p = ce ^ . If Q = qi/3, the lipids 
will follow a Boltzmann distribution under these assumptions, i.e., 



p = ce q '^. 



(2.11) 
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To solve for the constant c, we make use of the quantity 

T = [ pdS= [ ce~ qi ^ dS, (2.12) 



which is the total number of the charged lipids on the surface T. Since the lipids stay 
on the surface r, T is a conserved quantity, which gives a calculation for c as 

(2.13) 



e - « w dS 
r 



Plugging (2.13) back in for (2.11) gives a constraint Boltzmann distribution of lipids: 

p[V] = . (2.14) 

e- qi ^ dS 



This distribution depends on the position of T. 

2.2. Electrostatic potential energy. The following theorem establishes the 
form of the electrostatic potential energy (sometimes called the electrostatic free en- 
ergy) for the entire protein-membrane system. 

Theorem 2.1. The form of the electrostatic potential energy G is 

G[r ; ^]= ^ [£|V<£| 2 - + Xs B(<f>)]da+^(]n j e-^dS- In J dS^. (2.15) 

Proof. The result will be established by showing that the variational form of 
(2.15) is equivalent to the Poisson-Boltzmann Equation (2.5). 

The variational form of the Poisson-Boltzmann Equation (2.5), for an arbitrary 
test function ip : ft — > R is 

-V • {eV(j))xjj dX + I XsB' {cj>)4> dX = [ fipdX. (2.16) 
Jn Jn 

Splitting the first integral in the domains Q s , f2 m , and f2 p , and by the definition of 



-V • (s s V<p s )ip dX+ -V • {s m V<j) m )ip dX+ I -V ■ (epV^ip dX 
Jn m 

+ f B'((f>)ip dX = [ ftp dX. 



By the product rule for divergence, keeping in mind that is a vector and ip is a 
scalar, 

/ -V • (e s V0» dX+ I e s {\7<p s ) • (V^) dX + f -V ■ {e m V<j) m ip) dX 
Jn s Jn B Jn m 

+ [ e m {\7(f> m ) • (VV) dX + [ -V • {e p V(jfip) dX + [ e p (V0 p ) • (Vip) dX 



+ / B'(4>)TpdX= / fipdX 
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Combine the second, fourth, and sixth integrals into a single integral over f2 and use 
the divergence theorem on the first, third, and fifth integrals, keeping in mind that 
the multiple boundaries of the domains gives 



f e(V0) ■ (VV>) dX - [ s s ip(V(t> s )-ndS+ [ £ S V(V^ S ) ■ n dS 
Jn Jan Jr c 

+ [ £ s iiV<f) s )-ndS+ f e s ^(V(t) s )-ndS- f £ m ^(V0 m ) • n dS 

- / £ m V(V0 m ) -ndS- [ £ p V(V0 p )-nd5+ f B'{(j>)il> dX = [ /V dX. 
Jr c Jr p Jn s Jn 



But, the test function ip is supported compactly over f2, so the boundary integrals 
over dfl are zero. Combine the boundary integrals over T and obtain 



( / £(V</>) • (VV>) dX+ f B'{<P)^ dX - [ fip dX 
\Jn Jn s Jn 



[ I V(e s V0 s - e m V<j> m ) ■ndSj + (J i>(s s V<t> s - e m V<j> m ) ■ n dS 



+ (^J V(e s V0 s - e p V(j) p ) -nds\ = 0. 



Next, apply the boundary conditions to the last three integrals to obtain the varia- 
tional form of the Poisson-Boltzmann Equation, 



( f e{V4>) ■ (W) dX + [ B'(<j))ipdX- f fipdX) 
\Jq Jn B Jn ) 

p[T c ]qiip dS p\T e ]qiip dS = 0. 



(2.17) 



Next, we will see that the variation of the electrostatic potential energy is cquiv- 
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alcnt to (2.17): 



lim G[</> + ^]-G[0] = dG[t;<f>,il>] 



d 
dt 



7' 



dt 



Jr 



= / eV((t> + tip) • (W) +XsB'((/> + til})il> -/(</> + #) 



7 



/ 



-<nP(<t>+til>) 



e(V0) • (W) + X S 5'(«£)V - /V rffi + 



t=o 

T / e-^i-qiip) dS 



/ 



e(V0) • (W) + XsB'{4>)ip - fi> dfi - / p[T] qi iP dS. 



The computation utilizes the fact that ~^j\ x \ 



1 dx 

-. — -Ax - — ) for any vector a:. The result 
\x\ dt 



of the computation matches the form of (2.17), which finishes the proof. □ 

2.3. General lipid distributions. In the previous section, we assume the lipids 
distribution is governed by the surface electrodiffusion equation (2.6) so it follows the 
ideal Boltzmann distribution (2.14). The practical distribution of lipids is subject 
to various constraints such as finite sizes and entropy conditions, and will not follow 
this surface electrodiffusion equation [15, 14, 22]. The following theorem states that a 
more general distribution of lipids is allowed for our definition of electrostatic forces 
on membrane surfaces: 

Theorem 2.2. Suppose the lipids follow a positive distribution p of the form, 



p[r] = 



(2.18) 



qi j^l{4>) dS 



for some distribution 7(0). Then, the form of the electrostatic potential energy G is 



G[r ; , 



f4> + XsB{<t>) dQ + C* ( In J 7 (0) dS - ln J dS 



(2.19) 



Proof. Most of the details of this proof follow the logic of Theorem 2.1. The only 
differences arise from the second half of the electrostatic potential energy, which we 
shall call G2, i.e., 



G 2 [r ; <P] = C* (in J 7(0) dS - ln J dS 



(2.20) 
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Computing the variation of G2, 



lira 

t->o 



6*2 



til>]-G 2 [<i>] dG 2 {t-A^} 



dt 



t=o 



dn C*^p W dS-lnJ r dS 



t=o 



C* 



7 '(0 + #) dS 

r 

7 (0 + #) dS 



C*j'{4>)tp dS 



l(4>) dS 



t=o 



p[T]qiip dS, 



which matches exactly the corresponding second half of (2.17), finishing the proof. 
Note that for lipids with a Poisson-Boltzmann distribution, 7(0) = e~ qi ^^ and C* = 
T//3. Q 

3. Shape derivative of electrostatic potential energy. In this section, the 
main calculation of JrG[r] will be established via the method of shape derivatives. 
For the shape derivative calculation, a velocity function of the membrane movement 
will be defined. According to the Hadamard-Zolesio Structure Theorem of shape 
calculus [7], the variation of the electrostatic potential energy (2.15) with respect to 
the position of the smooth boundary T, that is, the shape derivative, must be given 
by the inner product of the force with the normal component of the velocity of the 
deformation. Thus, the shape derivative <5rG[r] will provide the force on a membrane 

r. 



3.1. Velocity of the surface. Define the velocity function V <E C°°(R 3 , 
the following dynamical system, 



V(x) Vi > 0, 
= X, 




3 )by 



(3.1) 



where X is the original position of the membrane and x is the transformed position. 
Assume that V is compactly supported near the bilayer membrane surfaces, i.e., 
V(x) = if dist(x, T "t Va ) > d for some d > where 



d < imin [dist (r e> flfi),dist (r c ,supp (/)),dist (r c ,supp (/))]. (3.2) 

This condition prevents the exterior membrane surface from stretching beyond the 
containment domain Q and also prevents the either membrane surface from overlap- 
ping the center of an atom contributing to the charge density / within the protein. 

The solution to (3.1) defines a diffeomorphism T t : M 3 — >• R 3 where T t (X) — 
x(t, X) maps the old coordinates X into the transformed coordinates x. By a Taylor 
expansion, we can approximate the map by 

T t (X)=x(t,X) 

= x{0, X) + td t x(0, X) + 0(t 2 ) 
= X + tV(x(0,X)) + 0(t 2 ) 

= X + tV(X) + 0(t 2 ), (3.3) 
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so that the map T t {X) agrees with the perturbation of the identity up to the leading 
term. 

The configuration under the transformation T t (X) will result in a new electro- 
static free energy, 



G[T t 



|V0| 2 - f4> + XsB((j>) dx + C* In / ~f(<p) dS -In dS , (3.4) 



where each of the functions are computed over the transformed regions fl t = T t (f2), 
r t = T t (T). Similar to the standard nonlinear Poisson-Boltzmann equation without 
surface charge distributions [16], there is a unique minimizer ip t & Hg(Cl) n X°°(f2) 
that minimizes (3.4) over Hg(£l). The minimum is 



G[T t }= min G[T t , 

<t>em(n) 



= G[T t ^ t \. 



(3.5) 



Following the ideas of Theorem 2.2, the same ipt is also the unique weak solution to 
the transformed boundary value problem of the Poisson-Boltzmann equation, 



0* 



[<t>] = 




on r Ct , r et , r Pt , 








s dn 


d<j> m rr , 
£ m dn p[T t \qi 


on r Ct , r et , 


d<j> s 
s dn 


dcf? 


on r Pt . 


£v ~dn~ 



(3.6) 



3.2. Transformation properties. The transformation T t (X) defined by (3.1) 
can act on volume and surface elements in tt. Some useful properties of the transfor- 
mation are outlined below, and are shown in details in [6]. 

3.2.1. Volume transformation properties. The following are assumed prop- 
erties of the transformation T t (X) defined by (3.1) on a volume element dX £ R 3 . 
These properties will be used in the computation of the shape derivative and their 
justifications arc found in [6], for example. 

(Tl) Let Xel 3 and t > 0. Let VT t {X) be the Jacobian matrix of T t at X defined 
by (VT 4 (X)). y = djT?(X), where T l t is the ith component of T t (i = 1,2,3). 
Let 



Jt{X) = detVT t (X). 
For each X, the function t H> Jt(X) is in C°° and at X, 

^ = J t (V-V)oT t . 



(3.7) 



(3.8) 



At t = 0, since no time has passed, VTo = / for any x and so Jo(X) = 1. The 
continuity of Jj at t = implies J t > for t > small enough. 
(T2) Define A(t) : O ->• R for t > small enough by 



A(t) = J t (VT t )- 1 (VT t )- T . 

At each point in tt, 

A'(t) = [((V • V) o T t ) - (VTt)" 1 ((VI/) o T t ) VT t 
- (VT)- 1 ((VV) o T t f (VT t )l A(t). 



(3.9) 



(3.10) 
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3.2.2. Surface transformation properties. In this section, we will compute 
a useful property of the transformation T t on a surface element dS €E K 2 . Define 

(3.11) 



J s (X,t) = (dctVT t (X))|VT t - J n(X)|. 



as in [5]. Note that at time t = 0, J s (X,0) = (det/)|/n | = |no(X)| = 1, since the 
length of the unit normal is always 1. Analogous to the transformation on a volume 
element, the differential surface element is transformed by ds = J s dS. The surface 
transformation property we wish to establish is the following: 
(SI) The derivative of J s in (3.11) is given by 



*£ = (j t (y.V)oT t )\\ 



— T 



dct VT t 
|VT t - T n| 



(VT ( 



— T„ 



dt^"^ 



(3.12) 



The above computation uses the product rule and (3.8). The remainder of this section 
will serve to complete the formulation of (3.12). 



First, the derivative with respect to t of (VT t T n) is given by 



Jt^ 



-n) = 



d(VT t " 



„ t dn 

VTr T — . 



(3.13) 



dt 1 dt 

There is an alternative expression to (3.13) that will prove easier for further calcu- 
lations. Let a = VT t _T n so that n = VT t T d. Then, computing the derivative of 
n, 



dn 
~dt 



d(VT t T d) _ d(VT t 2 



dt 



dt 



, da 
dt' 



Therefore, 



da 
"dT 



, VT ,_ T , /dn d(VT t T 
(VJ * j ldt dt 



Substituting back in the definition of a gives the alternative form of (3.13), 



-T 



dn d(VT t T ) 



dt 



(VT t " T )n 



(3.14) 



This equation is more advantageous than (3.13) because it only requires computing 
the derivative of VT t T , whereas (3.13) requires computing the derivative of VT t ~ T . 

The derivative of VT t T with respect to time t is next. It is important to keep 
in mind which variable, x or X, is the variable of differentiation when writing V(-). 
Sometimes the variable will be emphasized by writing V x (-) or Vx( - ) depending on 
the context. Write T t {X) = T(t, X) to emphasize T is a function of X and t (but not 
x). This will allow for the subscript notation of partial derivatives. Computing the 
derivative with respect to time of VxT(t, X) yields 



lv,T(t,X)=^.T i 



d x .Vi{X) = VxV{X). 
Now the formula V X V(X) = (V x V(x))(V x T(t, X)) will be verified: 



dx >^t 



V X V{X) 



( dV 1 
dX 1 

dV* 



dVi 
~dX> 



\dX 1 



dX 3 

dV 3 
dxj 



dxi 
\dx\ 



dVi 
dx2 



dx 3 

m 

dx 3 J 



( dx% 
dXx 

dx 3 
\dXi 



dxl 



dX 3 

dx 3 
dX 3 J 



(V x V(x))(V x T(t,X)). 
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This establishes 



dt 



V x T(t,X) = (V x V(x))(W x T(t,X)). 



(3.15) 



From (3.15) and the fact that the derivative of the transpose of a matrix is the 
transpose of the derivative, 



^V7f = ((W)(VT t )) T . 



(3.16) 



The equation (3.16) may be substituted in (3.14), which may in turn be substituted in 
to (3.12), but we choose to substitute everything at the end of the section for neatness. 
Next, we will establish a formula for the normal derivative % at time t = 0. Let 

' at 

u and v parametrize a surface in R 3 and define the surface by 

r(u, v) = (x(u, v),y(u, v), z(u, v)). (3.17) 



With this parameterization of the surface, the transformation T t acting on the surface 
can be defined as 

T t (x{u, v), y(u, v), z{u, v)) = (x(u, v) + a{u, w, t), y(u, v) + /3(u, v, i), z(u, v) + ~f(u, v, t)) 

= r{u, v ) + s(u, V, t), (3.18) 

where s(u,v,t) = (ct(u, v, t),{3(u, V, t),j(u, v, t)). Note that s depends in t but r is 
independent of t. Next, define the (initial) normal vector to the surface in the reference 
coordinates by 



i(X,0) 



\r u Xr v \' 



(3.19) 



where the subscripts denote partial derivatives. Then, at any time t > 0, the trans- 
formed normal vector is calculated by applying (3.18) to the initial normal vector 
(3.19), 



n{X,t) 



(r u + s u ) x (r v + s v ) 
\(r u + s u ) x (r v + s v )\ 



(3.20) 



For simplicity, let £ = (r u + s u ) x (r v + s v ). Next, the initial time rate of change of 
(3.20) will be computed: 



dn(X,t) 


.1 


VlCl) 




dt 


t=o dt 




t=0 



1 d( d ( 1 
]C\M +C Jt \](\ 



(3.21) 



t=o 



Define 



f = 7— — and 

ICI dt 



c- (- 

^dt vici 



The derivative — can be computed using the product rule for cross products. Notice 
dt 
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that r does not depend on t, the simplified form of P is 



P\t=0 = 4| ^(( r « + S «) X K + S v)) 



1 



-j7(r u + s u )j x (r v + s v ) + (r u + s u ) x (4:{r v + s v ) 



t=o 



1 



(r u + s u ) x (r v + s v )\ 



1 



\r u x r v \ 
The simplified form of Q is 

- t— ( — 

Qlt =° ~ Q Jt vie! 

= (r u x r v )- 



ds u 

— - — X 7*i> H - 7*7, X — - — 

dt dt 



— x [r v + s v ) + (r„ + s u ) x — 
dt dt 



ds, 



t=o 



= {{r u + s u ) x (r + «„)) |^p(C • 

N / rfS M dSy 



\r u x r„| 3 

This gives the desired result: 

dn(X,0) 1 /ds, 

dt 



(r u x r v ) • ( — x r v + r u x — 



, x r„ + r u x — - 
|r„ x r v \ \ dt dt 



+ (r„ x r«) 



-1 



k„ x r^| 3 



(r„ x r„) 



X 7\; "t - T*7i X 

" u dt 



To simplify, define 



ds v 

R= — xr v +r u x —. 
dt dt 



Then 



dn(X, 0) 
dt 



1 



\r u x r v \ 
1 

\r u x r v \ 
1 

\r u x r v \ 



R+(r u x r v ) 



-1 



\r u x r v \ 3 



((r u x r v ) ■ R) 



R- 



\r u x r v \ \\r u x r v \ 



%-R 



[R-n(X,0) (n(X,0)-R)] 



That is. 



dn{X, 0) 
dt 



1 



\r u x r v \ 



[R - n(X, 0) (n(X, 0) • R)} 



(3.22) 



(3.23) 



To finish the formula (3.12) requires substituting (3.16) and an expression lengthier 
than (3.22) into (3.14) and substituting the result into (3.12). The final expression for 

— at arbitrary t > would be exceedingly long, and the expression is unnecessary 
dt 

for the calculations of the shape derivative. Indeed, only the derivative of J s at t = 
will be used. To compute this quantity we notice that by (3.16) it follows 



-VT/ 

dt 1 



= ((W)(VT )) T = (W) T , 



(3.24) 



t=o 
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where the V operator acting on T t r (X) is with respect to X and the one acting on 
V(x) is with respect to x. Now, using (3.12), (3.14), (3.24), and (3.23) gives the result 



dt 



= (Jo(V-F)oT )|VT 




(VT " T n 



t=0j 



t=o 



= V • V- 
= V • V - 
= V • V 



<^>(§-^W> 



(dn(X,0) 



dt 



\r u x r v \ 



(R- n(n ■ R)) — n ■ (W) T n, 



where n = n(JT, 0) is the initial unit normal, and hence \n\ = 1. Notice that n(n ■ R) 
is the component of R in the direction of n. Then, R — n(n ■ R) is the component of 
R perpendicular to n. Thus, n ■ (R — n(n ■ R)) = 0. This leads to 



dJ s 
~~dt 



V • V - n ■ (W) T n. 



(3.25) 



Next, notice that for any matrix A and vector x, x ■ A T x = x ■ Ax. This reduces (3.25) 
to 



dJ s 
dt 



V • V - n ■ (W)n. 



(3.26) 



The right-hand side of (3.26) is exactly the surface divergence of V. This leads to the 
final expression of the initial time derivative of J s , 



dJ s 
~~dt 



(3.27) 



t=o 



3.3. Shape derivative calculation. The main theorem of the paper estab- 
lishes the shape derivative of the electrostatic free energy to the boundary, which 
corresponds to F n , the normal component of the dielectric boundary force. 



SrG[T] = -y |V^| 2 + ^|W> m | 2 - ejVC ■ n\ 2 + e m <y% ■ n)(V^ m ■ n) 
-B{4>v)-q lP [Y](Vr -n). 



(3.28) 



Let V £ C°°(K 3 , R 3 ) be a smooth map that vanishes outside a small neighborhood 
of the membrane surface T. That is, V(X) = if distpT, T) > d for some d > 
satisfying (3.2). Let the transformations T t (t > 0) be defined by (3.1). For t > 
the electrostatic free energy is given by (3.5), where the functional G[F(, •] is given in 
(3.4) and tpt is the weak solution to (3.6). For t = 0, the electrostatic free energy is 
given by (2.15) and ipo is the weak solution to (2.5). 

Theorem 3.1. Assume f G H 1 ^). Then the shape derivative of the electrostatic 
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free energy G[T] in the direction of V is given by 

5t, V G[T] =j(- yl v ^o| 2 + ^|VV'o m | 2 - e m |VV m ■ n\ 2 + £m (V^ • n)(V^ • n) 
- B(i> ) - qip[T](V% ■ n) ) (V ■ n) dS. (3.29) 



Proof. The proof is divided into four steps. 

i. First, the energy functional will be computed in the transformed coordinates 
through a new function z(t,<fi). A change of variables will bring z back to the 
reference coordinates, and then z will be differentiated with respect to time. 

ii. Second, the difference quotient corresponding to the shape derivative will be 
squeezed between two realizations of d t z. 

iii. Third, the inequality will be passed to the limit as t — > and it will be shown 
that the two realizations of d t z are identical in the limit, and hence equal to the 
shape derivative. 

iv. Fourth and finally, the result will be simplified to match the final form. 

The computations to determine the shape derivative of (2.15) will be completed 
in two calculations. Splitting (2.15) into two components, 



Gi[T;<P}= f ^|V0| 2 -/0 + Xs S(0)l dn (3.30) 



and 



G 2 [r;<t>] = j (in J e- qi ^ dS - In J dS^j (3.31) 

where G[T; <fi] = Gi[T; <fi] + <fi]. The details of the computations for Gi[T; <fi] are 

omitted as they appear in [16] with slight sign modifications. 
Step 1. Let t > 0. 

Since each <f> <S Hg($T) corresponds uniquely to o T t -1 g Hg(Q), by (3.5), 

^1 = ^ qr^oTr 1 ]. 

The purpose of this is to take the transformed coordinates x = T t (X) back into the 
original coordinates, for which information is known. Let <f> <E iJ 1 (57) n L°°(^X) and 
t > and denote 

z{t,4>) = G[T t , ( f>oT t - 1 }. (3.32) 

The function z will be split into two components, zi, corresponding to G\, and Z2 
corresponding to G2. 

z 1 (t,4>) = -G 1 [r t ,<t>aT t - 1 ], (3.33) 
z 2 (t,4>) = G 2 [T u 4>oT-\ (3.34) 

so that z(t,(f>) = —zi(t,(f>) + Z2(t,4>). The reason for the sign change on (3.33) is to 
match the analysis in [16], which uses energy maximization, and here we consider 
energy minimization. Since the details of —zi(t,(f>) appear exactly in [16], we will 
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focus on only z 2 (*, <f>) here. Using (3.34), (3.31), and a transformation of coordinates 
x = T t (X), 

Z2(t,4>) =G 2 [T t , ( j ) oT t - 1 ] 

= C* (hi J j((f> o T t _1 )(x) dS(x) - In J dS(x)^j 

= C* (\nj^ j{cj>(X))J a (X,t) dS(X) -In J J s (X,t) dS(X)^j , (3.35) 

where the formula for J s is given by (3.11). Note that <f>(X) as a function of X does 
not depend on t. Differentiating with respect to t, 



d t z 2 (t,4>) = c* 



\ 



f ^ dJ s {X,t) 
7W j t dS{X) 

j(4>)j s (x,t) ds(x) 



dJ a (X,t) 
dt 



dS(X)\ 



J s (X,t) dS(X) 



(3.36) 



Then, the full form of dtz(t,4>) is given by combining the calculation in [16] and 
(3.36), 



d t z(t,4>) 



\A\t)V<t> • - ((V ■ (/V)) o T t )4>J t + x»B(0)((V • V) o T t )J t 



dX 



c* 



\ 



7 (0^M dS (x) 
j{4>)j s (x,t) ds(x) 



r 



dJ s (X,t) 



dS(X)\ 



J a (X,t) dS(X) 



(3.37) 



where A'(i) is given by (3.10). 

Step 2. Let i g (0,r]. Since t/>t € H L°°(Q) and Vo € ffi(fi) n L°°(Q) 

minimize G[rt, •] and G[r, •], respectively, over Hg{Vl) we have 



By (3.38), 



G[T u 4^oT^] > G[r t ,^ t } = G[r t ] 
g\t,1h] > G[r,Vo] = G[r] 
G[r, ^ o r t ] > G[r, Vo] = G[r] 



G[r f ,v.oor t - 1 ]-G[r,v.o] > G[r f ] - G[T] 



(3.38) 
(3.39) 
(3.40) 



By (3.40), 



G[r t ] - G[r] G[rvA t ] - G[r, ^ o r t ] 



t ~ t 

Putting the two inequalities together with the definition in (3.32) gives 
z(t, - z(0,^o) > G[r f ] - G[T] > z(t, ^ o T t ) - z(0, ^ o T t ) 



t 



t 



t 



(3.41) 
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Notice that the far left expression of the inequality is the secant line of z(-,ipo) 
from to t < t, and the far right expression is the secant line of z(-,ip t o T t ) from 
to t < t. Since z is differcntiablc on t > 0, we can apply the Mean Value Theorem to 
each expression. That is, there exists £(t),r}(t) € [0,t] for each t S (0, r] such that 

dtzM), i>o) > G[Tt] 7 G[T] > d t z(r}(t), ^ ° T t ) Vte(0,<r]. (3.42) 



Step 3. Claim: 



lim9tz(^(t),V'o) =9 t z(0,ip ), 



t->o 



Bm^(77(t),^toT t ) = d t z(0,i) Q ). 



t->o 



(3.43) 
(3.44) 



As in Step 1, the proofs for (3.43) and (3.44) will be split into the corresponding 
Z\ and Z2 components, 



\und t zi(£_(t),ifj ) = dtzi(0,vlx)), 
limd t zi(r](t),ipt ° T t ) = d t zi(0,ip ), 
\ymdtz 2 {£,{t),ip ) = 9 t z 2 (O,-0o), 
\imd t Z2{r)(t),^t °T t ) = d t z 2 {0, ip ). 



(3.45) 
(3.46) 
(3.47) 
(3.48) 



The justification of (3.45) and (3.46) are found in [16], so (3.48) will be proven 
here, noting that (3.47) is similar. Let rj(t) € [0, t] and do consider passing dtZi{r}{t), ip t o 
T t ) to the limit as t -> 0, 



d t z 2 (n(t),i> t oT t )=C* 



dt 



y(i> t oT t )J a (X,r)(t)) dS(X) 



'> dS(X)\ 
) 



J a (X,T)(t)) dS(X) 



c* 



( I 7(V>o) 

J To 



dJ s (X, 0) 
dt 



dS(X) 



dt 



-y(i/> o )J a (X,0) dS(X) 



'> dS(X)\ 
I 



J B (X,0) dS(X) 



By passing (3.42) to the limit as t — > and using (3.43) and (3.44) through the 
computation (3.37) shows that 

S r ,vG[r} = Jim G[Tt] ~ G[T] = d t z(0, Vo) 



-a'(o)Wo ■ v^ - (v • (/y))Vo + xsBWo)(y ■ v) 



dX 



c* 



, dJ s (X,0) ^ ofv ^ 

TW) "77 <^(X) 

r rf< 



7 (^o)J s (^,0) dS(X) 



dJ a (X,0) 
dt 



dS(X)\ 



J s (X,0) dS(X) 



(3.49) 
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Step 4. The final step is to simplify the form of (3.49) so that it matches the 
right-hand side of (3.29). The simplifications will begin with Z\ and then in z%. The 
combined results will be the shape derivative. In [16], c\zi(0, is simplified to 

d tZl (0,iP )= ( £ -±\vr \ 2 (V-n) dS- f £ -^\V^\ 2 {V-n) dS 



- J e s (Vr ■ n)(V ■ Vr ) dS + J e m (V^ 1 ■ n)(V ■ V^) dS 

+ [ B(ipo)(V ■ n) dS. (3.50) 



The tangential components of an d ip™ are equal on any interface T, but the normal 
components of course are not, so any difference occurs only in the direction of n, 

V(V>o - WT) = (Wo • n - V^J 1 • n)n on T. (3.51) 

Now apply the new interface condition: 

e s ■ n = e m VV-o™ ■ n - qip[T] , (3.52) 

where p[T] follows (2.18). There are two equivalent approaches from here. Consider 
the third and fourth terms of (3.50). Substituting (3.52) into the third term and using 
(3.51) gives 



e,(WS -n)(V- V%) dS + J e m (V^" • n)(V ■ V^ n )dS 

- f (e m V$P • n - qip[T])(V ■ W%) dS + J ^(Vff • n)(V ■ V?/> m ) dS 

■(e ro V^ • n)V ■ V(^ - % n ) + qip[T](V ■ V%) dS 

-(e m W^ ■ n)V ■ (V^g • n - • n)n + qi P [T}(V ■ V^o) dS 

m • n| 2 (V • n) - e m (V^ ■ nXWS" • n)(V ■ n) + qip[T](V ■ V^) <*& 

Substituting this term back into (3.50) and combining terms will lead to the final 
result for dtZ\, 



d t zi(0,ijj ) 



1 1 iVVof (^ ■ n) dS - jf ^ W| 2 (^ • n) dS 

m |V$? 1 • n| 2 (V • n) - £ m (VVo • n)(V$? • • n) + qip[T](V ■ V^) ^ 
+ J B{i/fo)(V-n) dS 

= J [(y |V^| 2 - ^|V^| 2 + UVVo m ' ^l 2 " ' »)W ' «) 

+ B(V>o)) (y-n) + ^[r](y-V^)] ^ (3.53) 

An alternative form of (3.53) can be obtained by substituting (3.52) and keeping 
(e s X7ipQ ■ n) instead. 
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Next, the shape derivative from d t z 2 will be simplified. Evaluating (3.36) at t = 
and 4> = "00 gives, 

dJ ^ dS(X)\ 



d t z 2 {QM) = C* 



7(V>o) j t dS{X) 



dt 



\ 



7 (V-o)J s (^,0) dS(X) 



J s (X,0) dS(X) 



The time derivative of J s at t = is the surface divergence of the velocity, as calculated 
in (3.27). Recall that J S (X, 0) = 1. Using these two facts, the above equation is 
equivalent to 



I I 7(V'o)(V s • V) dS(X) [ (\7 s -V)dS(X)\ 



d t z 2 (0^ ) = C* 



\ 



7 (Vo) dS(X) 



dS{X) 



(3.54) 



/ 



Next it will be shown that the surface divergence of any continuous function on 
an enclosed surface is 0. Let F be an arbitrary function and consider the following 
integral on the arbitrary enclosed surface Fq, 



/ (X-F)dS 2 (X). 
Jr 



(3.55) 



The notation is suggestive of application to the problem in context, but is for now 
treated as arbitrary for the purpose of generalizing. The notation Si(X) indicates a 
surface in W. Applying the divergence theorem, 



/ (V ■ F) dS 2 {X) = f F-n 1 dS 1 {X)+ I F ■ n 2 dS l 



(X). 



(3.56) 



Since Tq is an enclosed surface, the outward normal directions to the boundary 8Tq 
are opposite in sign. Therefore, the two integrals in (3.56) cancel, 



/ (V s -F) dS 2 (X) = 0. 



(3.57) 



By applying (3.57) to (3.54), the second term is zero, as the velocity V is contin- 
uous. Applying the product rule to the numerator of the remaining term, 



d t z 2 {0,ip ) = C* 



■ (V 7 (ip )) - V ■ (V s <y(i> )) dS(X)\ 



r 



V 



7 (Vo) dS(X) 



(3.58) 



r 



/ 



By applying again (3.57) to (3.58) to remove the surface divergence of the continuous 
function V~/(ipo), and by expanding the surface gradient in the remaining term, we 
have 



d t z 2 (0,Vo) =C* 



-V ■ ( 7 '(^ )V s Vo) dS(X) 



7 (Vo) dS(X) 
q lP [T](V -X s ^o) dS{X). 



(3.59) 
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where the last step uses the definition of p[T] in (2.18) at = ipQ. 

Finally, combine (3.53) and (3.59) for the full form of the shape derivative, 



S r , v G[T] = -d tZl + d t z 2 

"(-y IVV^oI 2 + ^ W| 2 - e m \V^ ■ n\ 2 + e m (V^ • n)(V^ m • n) 



BWo)) (V-n)-q lP [T](V-Vr ) 



dS 



qi P [T}(V -V s ^ )dS. 



Notice that ipQ = V>o' on T by the first boundary condition of (2.5), and hence the 
last two terms reduce: 



-qip[T](V -W%)dS{X)+ / qip[T](V- V s Vo) dS 



qip[T]V ■ (VWo - VV»o) dS 

q lP [T]V ■ ([Wo - (VV^o ' n)n] - V^) dS 



- qi p[T}(Vr -n)(V-n) dS. 



Using this result, the shape derivative is simplified to the final result, 



B(^0) - gip[T](WS • n) (F-n) dS. 



(3.60) 



by 



The alternative form of <$r,vG[r] based on the alternative form of (3.53) is given 



4 al ^G[r] = 



B(^o) ~ qip[T](y% n ■ n) ) (V ■ n) dS. 



(3.61) 



Notice that the two forms are equivalent, since subtracting (3.61) from (3.60) using 
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the second interface condition in (2.5) gives 
8 rtV G[T}-S^G[T} = 

= y(- e TO |V$P • n\ 2 + e m (Vr ■ ™) W • n) - «p[T](V^ • n)) 

- (e s |V^ • n| 2 - e ,(V^ ' «)(VC ■ n) - ffi/>[T](VC ■ n)) j (V ■ n) 

= ^ - EmlVVS* • n\ 2 + e m (V$5 ■ n)(V% n ■ n) + (e s (V^ ■ n) - e m (VC ■ n)) (v^ • n) 

- e.lV^S ' "I 2 + £,(V^o ' n)(W ■ n) + (e m (V!T • n) - e.(V^S • n)) (v^ ■ n) J (V • n) 
= 0. 

□ 

4. Equivalence to the Maxwell Stress Tensor. Define E = —Wtp. Then the 
Maxwell Stress Tensor is defined by 

T = eE®E- £ -\E\ 2 I-XsB{4>)I. (4.1) 

We shall verify that the force F n given by (3.29) can also be expressed as the jump 
in the Maxwell Stress Tensor [1], 

F n = n ■ T + n — n ■ T~n 

= (e s \Vr ■ n\ 2 - £ f\Vr\ 2 - XsBm) - (e m \V^ m ■ n\ 2 - £ -f\^ m \ 2 ) 

= -y IVV1 2 + ^Wf + e.W ■ - e m \Vr i ■ n\ 2 - B{^) 

= -y IVV1 2 + ^IW>"T + e s (S/r ■ n)(Wr ■ n) - e m (V#™ • n)(V7/> s • n) 

+ e m (\/r ■ n)(V^ m • n) - e m (V^ m ■ n)(Vf ■ n) - B(^) 
= -y IVV1 2 + y^VVH 2 - e m |VV> m ■ n| 2 + e m (VV> s ■ n)(W m • n) 

- B{i>) + (e s {Vr ■ n) - e m (VV> m • n)) (VV> S ■ n) 

= -y|VV s | 2 + ^Wf - e m |V^ m ■ n| 2 + e m (VV> s ■ n){V^ m ■ n) 

-Bdj)-q lP [T}(Vr -n). (4.2) 

This matches (3.28) exactly, suggesting that our definition of the electrostatic forces 
on charged dielectric interfaces is physical. 

Acknowledgments. The authors thank Benzhuo Lu for many helpful discus- 
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